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■ We extend our exact reformulation of the bosonic many-body problem in terms of a stochastic 

' Hartree ansatz to a stochastic Gutzwiller ansatz for the Bose-Hubbard model. The use of this ansatz 

^s^j , makes the corresponding Monte Carlo scheme more efficient for strongly correlated bosonic phases 

like the Mott insulator phase or the Tonks phase. We present a first numerical application of this 
stochastic method to a system of impenetrable bosons on a f D lattice showing the transition from 
the discrete Tonks gas to the Mott phase as the chemical potential is increased. 
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I. INTRODUCTION 



The Bose-Hubbard model of lattice bosons with on-site interactions Q is a widely used theoretical model for the 
description of a number of different physical systems, ranging from 4 He adsorbed in porous media Q to granular 
' superconductors and Josephson junction arrays [51 Irt m IsL 1^. ITfij . 

In the last years, the Bose-Hubbard model has attracted a renewed interest since its experimental realisation with 
ultracold atoms loaded in the periodic potential of an optical lattice 0, ^| . Following the suggestion of 01 1 this 
1 system has led to the observation of the superfluid-Mott insulator quantum phase transition and of quantum 

revivals of the matter field 0] . In the present paper we propose a new tool for the study of the Bose-Hubbard model 
£Nl ■ in the form of a Quantum Monte Carlo scheme based on the stochastic evolution of a Gutzwiller ansatz. This new 
scheme generalizes our previously proposed schemes based on a stochastic Hartree ansatz ^) . 

Intuitively, it is clear that the Hartree ansatz is not the optimal one for the study of a Mott insulator state: the 
I expansion of a Mott-insulator state on Hartree states involves an average over all the possible relative phases between 
£f) . adjacent sites, so that a large number of realizations is required for an accurate sampling of the state. Furthermore, 
if one rises the on-site interaction strength to infinity so to realize the so-called impenetrable boson model, the rate 
of increase of the statistical error with time, being proportional to the strength of the atom-atom interactions, will 
1 diverge. 

On the other hand, approximate analytical mean-field calculations based on a Gutzwiller-type wavefunction have 
given predictions which are in reasonable agreement with the exact numerical calculations jl8j |. This suggests that 
a suitable stochastic generalization of the Gutzwiller mean-field equations could lead to an efficient quantum Monte 
Carlo scheme. The main difference with respect to the reformulation in terms of a stochastic Hartree ansatz consists 
in the exchanged roles of the kinetic and the interaction energies: in the Gutzwiller case, interactions can be fully 
taken into account by the deterministic evolution, while the noise is necessary to keep track of the correlations 
between different sites which are induced by the kinetic energy term in the Hamiltonian. While the stochastic Hartree 
reformulation was well suited for the study of weakly interacting Bose gases |l9j . the present Gutzwiller scheme is 
expected to be more suited to the study of systems of strongly interacting bosons, such as the Mott phase or the gas 
of impenetrable bosons. 

In sec|n] we review the Bose-Hubbard Hamiltonian and we introduce the Gutzwiller ansatz. In sec lllll we show 
how it is possible to find a suitable stochastic evolution of the Gutzwiller wave function such that the stochastic 
evolution of the ansatz is equal, in the average, to the exact evolution given by the Bose-Hubbard Hamiltonian 
In sec lIVI we discuss the evolution of the statistical error of this stochastic reformulation of the many-body problem 
and, in particular, we show that the error remains finite at any time of the evolution. This guarantees that the 
stochastic Gutzwiller reformulation can be actually used for numerical Monte Carlo calculations. In secE] we extend 
the reformulation to the case of an imaginary-time evolution so as to obtain the thermal equilibrium state of the 
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Bose-Hubbard system at a given temperature and chemical potential. A first application to the case of impenetrable 
bosons on a ID lattice is presented, which shows the transition from a discrete Tonks gas to the Mott phase as the 
chemical potential is increased: the results are in good agreement with analytical calculations performed by means of 
fermionization techniques. For the sake of completeness, a brief discussion of these techniques is given in the appendix. 
In sec lVII we generalize the results to the case of several spin components. Conclusions are drawn in sec lVIII 



II. THE MODEL HAMILTONIAN AND THE GUTZWILLER ANSATZ 



In its simplest form, the Bose-Hubbard Hamiltonian can be written as: 

H = Hu+Hj = jY. * f (r)^r)*(r)*(r) + ^ £ <5« *t(r) *(r') , (1) 

r r.r' 

in which the spatial coordinate r runs on a D-dimensional cubic lattice of Af = Y\f = i A/i points with periodic boundary 
conditions. With li and Li = Nik, we shall denote respectively the lattice spacing and the total lattice size along 
the i-th dimension. The on-site interaction is taken into account by the first term TLjj, which is proportional to the 
charging energy U . The hopping from one site to a neighboring one is described by the other term TLj. J quantifies 
the amplitude of the hopping processes. The function 8^, is equal to one if r and r' are neighboring sites, otherwise 
it vanishes. The annihilation operators ^(r) satisfy the bosonic commutation rules [^(r), ^(r')] = S r r >. 
The generic Gutzwiller ansatz can be written in terms of the complex wavefunctions <p(r, n) as: 

OO -j- 4- / \ ^ 

|G:0)-nE^ r ' n )— 7=- l°>- (2) 

r n=0 v 

This kind of states has been extensively used within mean-field theory to obtain approximate analytical results for 
strongly interacting bosons on a lattice: since the wave function has a factorized form, no correlations between different 
sites are included and the effect of Tij can not be treated exactly. This is complementary to what happens in the 
case of a Hartree ansatz where the mean field theory can deal exactly with the hopping Hamiltonian, but not with 
the interaction one. 

The generality of the Gutzwiller ansatz is however large enough to include Glauber coherent states as well as Fock 
states. A coherent state of macroscopic wavefunction 7(r): 



oo 



|coh : 7 ) = e-'HI 2 / 2 ]T -L(V 7 (r)*t (r) ) n | ), (3) 

where ||7|| 2 = ^ r l7( r )| 2 ' corresponds to a Gutzwiller state with: 

0(r,n)= e -W r )l 2 / 2 ^. (4) 



On the other hand, a Fock state in which each site has a well defined occupation no(r) corresponds to a Gutzwiller 
wavefunction of the form: 

c/)(r,n) = <5„,„ ( r ). (5) 

The Gutzwiller ansatz is therefore much more flexible than the coherent state ansatz often used in quantum optics. 
Its normalization is given by: 

oo 

(G:0 s |G:0 d )=nE^( r ' n )^( r ^)- ( 6 ) 

r n— 



HI. REAL TIME EVOLUTION 



As the set of all Gutzwiller states forms an overcomplete family, one can always assume that the initial state vector 
IV'(O)) is a coherent superposition of Gutzwiller states. As discussed in 01 f° r the case of a Hartree ansatz, one can 
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reabsorb the phases of the complex coefficients of this superposition by a redefinition of the Gutzwiller wavefunctions, 
so that the initial state vector can be written as the average of the Gutzwiller ansatz over some probability distribution 
for <i>: 



\G 



(7) 



where we represent this average by angular brackets (■••)• In the present section, we shall show that a suitable 
stochastic evolution can be found for the Gutzwiller wavefunctions <p(r,n), such that a relation like (JJJ) holds at any 
time t > 0, in the sense that the exact state vector at time t is given by the average over the ensemble of Gutzwiller 
states resulting from all possible stochastic evolutions of the initial ensemble during the time t: 



\m) = (\G: m 



(8) 



In the hopping Hamiltonian Tij, we first explicitely isolate the mean- field contribution by rewriting: 



Hj= 3 - J2 5 ?r'* f (r)^(r') - J - [(* t (r)-^*(r))(*(r')-Vi(r')) + *(r)^(r') + * t (r)V;(r')-^(r)^(r') 



in terms of the C-number mean-field ^(r)- defined as: 

<G:0|*(r)|G:^) _ (0(r)| * |0(r)> 



(9) 
(10) 

(11) 



^ (r) (G : t\G : </>) (#r)Wr)> 
where the ket |0(r)) describes the state of the site r: 

oo 

|#r)) =J2<t>M\n). 
n=0 

This corresponds to the general strategy, already implemented for the Hartree-Fock ansatz for bosons 01 an( l f° r 
fermions [20|, to have a deterministic evolution identical to the mean- field theory so that the stochastic part includes 
only deviations from mean-field. Mathematically, as we shall see, this makes the noise term strictly orthogonal to the 
state \G : <fi), which is expected to minimize the growth rate of the statistical error. Note that only the first term 
of the square bracket JUJ can induce correlations between different sites, and it is this term which will be taken into 
account with the noise term of the evolution. 

The action of the total Hamiltonian Q on a Gutzwiller state of the form J2J) can then be written as: 



±H\G^)dt = Y,[ E X( 



r; n, n m r, n 



&t( r )n 



r n,n'=0 



[nxx^ 



#t(u) n 



|0) + 



Jdt 
Yi 



, , \ #t( r )«^t( r ')«' 



in — l \ / \ / \n\rv\ 

n,n'—0 

[ n f>(u,»)* t(n) " 



u/r,r' n— 



where the matrices %(r; n, n') have vanishing entries for \n — n'\ > 1, 

dt\Un(n-l) J 



x(r; n, n) 



27V 



10), (12) 



(13) 



X (r;n,n+1) = ^ Vn + 1 ^ *S^V) 



2^ 



(14) 



and 



X(r; re, ra - 1) = ^ ^ <^ r W(r'^ 



(15) 
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so to include all the on-site terms coming from both the interaction and the mean-field part of the hopping Hamiltonian 
©. The factor 1/7V appearing in the last term in H13[) comes from the normalization 10 of the Gutzwiller state. 

We now consider a generic stochastic evolution for </>(r, n): 



dcj>{r,n) = d4> d (r,n) + d4> s {r,n) 



(16) 



where d<f>d(r, n) is a deterministic drift term proportional to dt and ckf) s (r, n) is a stochastic noise in the Ito sense. In 
particular, the mean value of d4> s is zero and its variance proportional to dt. The average evolution of the Gutzwiller 
ansatz during dt can be obtained by means of an expansion in powers of dt; retaining the terms up to order dt, one 
is left with: 



G 



G ^> + E [ E d Mr, n)^} \ J] £ n)^&£] |0> 
_ Vn\ J kti„ n vn! J 



n=0 



u/r n— 



|E[f < %(,,»)^. B -)* , ">y] [ n £«u,n)*^]|0> (17) 



u^r,r' n— 



The stochastic evolution (|17|l can be made to coincide with the exact one given by (|12[1 simply by imposing that the 
deterministic evolution d(f>d is equal to: 



dt 



b d (r, n) = E X(r; n, n')(f>(r, ri) = — 



Un(n - 1) 



[r,n 



E <© W) n - 1) + ^*(r') V£ + T#r, n + 1)) - ^ ( E *£ ^>Mu')) #r, n) 



(18) 



and that the correlation function of the noise term satisfies: 

d^(r, n) d0 fl (r',n') = ^r^J [(Vn^(r,n - 1) - $*(r) 0(r,n))(V^+10(r' ) n / + 1) - ^(r') 0(r',n')) + 



+ (v / rT+10(r,n + l)-V'(r)<?!)(r,n))(V^ 7 ^(r / ,n'-- 1) - V*(r') 00', n')) • (19) 

The symmetrized form of (|19|) with respect to (r, n) and (r', n') is imposed by the general property of a noise correlation 
matrix being symmetric. Apart from an additional term giving a mere phase rotation, the deterministic evolution 
(|18H coincides with the usual evolution considered in mean-field theory. It is easy to verify by direct substitution that 
a noise with the required correlation function (|19f) can be obtained as: 



d(f> a {r,n) = E 



i w(k) 
2Af 



1/2 



dB k e lk ' : 



Vn0(r,n - 1) - ip*(r) <j>{r, n) ) + 



+dBt e- lk r 



Vn+T 4>{r, n + 1) — ip(r) 0(r, n) 



with c?_Bk being independent and zero- mean complex random variables such that: 



dB k dB k , = 



dB k dB£, = 2dt5 K v 



(20) 



(21) 
(22) 



For the lattice here considered, the quasi-momentum k is a .D-component vector whose ith component is equal to 
ki = 2Tmi/Li, rii being an integer number < < Afi- w(k) is the particle dispersion in the lattice in the absence 
of interactions: 



fiu;(k) = J^cos(fcjij). 



(23) 



In any practical simulation, a cut-off -/V max has to be set in the sums over the site occupation numbers, which physically 
corresponds to projecting the dynamics on the subspace in which all sites have an occupation number smaller or equal 
to iV max . The error induced by imposing this cut-off is negligible provided a sufficiently large value of -/V max is adopted. 
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IV. EVOLUTION OF THE STATISTICAL ERROR 



As we have discussed in detail in 'JJ% , formally imposing the correct value of the averages is not in general sufficient for 
a practical application: it is essential to prove the finiteness of the statistical error, which guarantees the convergence 
of the Monte Carlo simulation to the exact result in the limit of an infinite number of independent realizations. 

As usual, the variance of the statistical error on the state vector of the system is defined in terms of the stochastically 
evolving Gutzwiller ansatz \G : 4>{t)) and the exact state vector \i(>(t)) as: 

E(t) = (\\\G:cj ) (t))-\m)f)- (24) 

Using the definition (jHJ of the stochastic reformulation, one can write the evolution of \G : (j>(t)) during an interval dt 
as: 

d\G : 4>(t)) = ^U\G : </>$)) + d\G : </>(t)) „ (25) 

where the Ito noise term d\G : <f>(t)) depends linearly on the stochastic term d(j) s (r,n) of the evolution of the 
Gutzwiller wave function. Our specific choice of noise terms H20fl not only leads to the orthogonality of this Ito 
term to the Gutzwiller ansatz: 

(G:ct>\ d|G =0 (26) 

but, even more strongly, to a purely orthogonal stochastic variation of the state of each site r: 

oo 

^>*(r,n)d<Mr,n) = 0. (27) 

Using the fact that the square modulus of the exact state vector HIV'W)!! 2 is constant and equal to 1 during the 
evolution, the variance of the statistical error can be rewritten as: 



E(t) = (A(t)) 1 (28) 

with: 

A(*)= |||G:0(i))|| 2 . (29) 

Inserting the explicit form of the stochastic noise term (|2U|I and taking advantage of l|27|) , one obtains an upper bound 
to the variation of dA in the interval dt: 

dA = \\ d\G : 4>) s\\ =A(f)2^ (0( r )|0(r)) ^~t J U + 1 ) A (^ ( 30 ) 

and therefore on A(t) itself: 

A(t) < A(0) e ^Dj( NlDax +i)t/n_ (31) 

This shows that the stochastic trajectories for the Gutzwiller ansatz cannot escape to infinity in a finite time. By 
ensemble averaging, one immediately finds an upper bound to the variance E(t) of the statistical error l|28|l: 

E{t) < [E(0) + 1] e MDJ ^ N ^ + ^ t ' h - 1. (32) 

This inequality guarantees that the statistical error of the stochastic reformulation using Gutzwiller states remains 
finite at all evolution times. 

This result is to be compared with the one we have obtained in [l5| for the Fock simple scheme: 

E(t) < [E{0) + l]e NUt/h -l, (33) 

where TV is the number of atoms. It is apparent how for TV sa J\f and -/V max w 1, the present scheme is more suited 
to simulate a system in the Mott-insulator phase U J, for which the interaction energy dominates over the kinetic 
energy. 
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V. IMAGINARY TIME EVOLUTION 



The previous sections have been devoted to the discussion of a stochastic reformulation of the real-time evolution of 
the Bose-Hubbard model. On the other hand, the density matrix p — exp(— H/ksT) describing the thermodynamical 
equilibrium state of a system at a given temperature T can be obtained by means of the imaginary-time evolution: 



dp 

~d7 



1 



{H,p(r)} 



(34) 



during an imaginary-time interval r = — > (3 — 1/kgT starting from the identity p(0) = 1 at r = 0. The expectation 
value of any operator O can be then obtained from p{(3) as: 



Tr[p6] 
~ 7 MpJ 



(35) 



As we have done in |16| for the Fock simple scheme, the stochastic reformulation in terms of a Gutzwiller ansatz 
originally developed for the real-time evolution can be easily extended to the imaginary-time evolution. 
It follows from the completeness of coherent states that the identity matrix can be expanded as: 

(36) 



(37) 



P(°) = 1 = ^77 J Vr y |coh:7)(coh: 7 | 
in terms of coherent states of the form J2J. The integration element T)j is defined as usual as: 

= Yl dRe[j{r)] dlm[ 7 (r)] . 



According to l|36l) . the imaginary-time evolution therefore starts from initially identical Gutzwiller wavefunctions 
^>i j 2(r,n;0) of the coherent state form J3J with a randomly chosen 7 (r) and proceeds as described by a pair of 
stochastic differential equations for 4> a =i t 2'- 



b a (r, n) = d4> d , a (r, n) + d4> s>ci (r, n) . 



(38) 



Apart from the obvious substitution % — > — the imaginary-time stochastic equations have exactly the same form 
as l|18ll and H20|) . As the Gutzwiller ansatz does not have in general a well defined number of particles, it is more 
convenient to use the Grand-Canonical ensemble with a chemical potential p, so that H = H — pN, N being the 
operator giving the total number of particles. The deterministic part of the evolution then reads: 



dr r fU n(n — 1) 



pn)(f) a (r, n) 



J 



2 J2 S rr> (^a(OVn0a(r,n-l)+^*(r')v^ + T^ a (r,n+l)) 

r' 

- ^(E6*;(*m) ^( r '«)] ( 39 ) 



and the stochastic part: 



dB k e lk r (VEcf>{r, n - 1) - i>* a (r) Q (r, n) ) + 

+dBle~ lUr Ufn + l(j) a {r,n + l) - ip a (r) <f> a (r, n) 
with dBk being independent and zero- mean complex random variables such that: 



dB k dB k , 







dBk dB^, = 2 dr <5k,k' • 



(40) 



(41) 
(42) 



In the practical simulation presented below, the dB^ have a Gaussian distribution. For the left and right Gutzwiller 
ansatz a — 1,2, ^> a (r) is defined as the mean-field corresponding to the wavefunction 4> a (r, n): 



V'c(r) 



(G:0 a |tt(r)|G:&,) 
(G : (j> a \G : <t> a ) 



(43) 
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The finiteness of the statistical error at any time of the imaginary-time evolution can be proven by means of the same 
arguments invoked above for the case of a real-time evolution. The evolution of the state vector during dr is given 
by: 

dr 

d\G : 0(r)) = -yff|G : 0(f)) + d\G : <j>(r)) s . (44) 

As the spectrum of the Bose-Hubbard Hamiltonian H is bounded from below in the presence of the cut-off at N max , 
an upper bound can be set to the growth rate of || \G : 4>(t)) || 2 , which guarantees that the statistical error remains 
finite at any time. 




' — — 1 1 1 1 1 1 1 1 1 1 

-3-2-10 1 2 3 



M/|J| 



0.4 




_Q -| I i i i I i i i I i i i I i i i I i i i I i i i I i i i I i i i 

'-16 -12 -8 -4 4 8 12 16 

X/l 



FIG. 1: For a ID lattice with 32 sites and the impenetrable boson regime U = +00, (a) mean filling factor n as a function of the 
chemical potential, (b) for fi = — |J|/2, first order correlation function g^(x) = {^(^^(O)) (open squares) and second order 
correlation function g^ 2 \x) = {^\x)^\Q)^> (0)^> (x)) (black disks) as function of the coordinate x. In (a) and (b), the solid 
lines are the analytical results obtained by means of fermionization tecniques as detailed in the Appendix. In both panels, the 
temperature is related to the (here negative) hopping constant by | J|/fcsT = 9/4. In (b) the difference between the chemical 
potential and the single particle ground state energy J is (9/8) &bT so that the fermionic analog of the Bose system is on the 
border of being degenerate. The position x is in units of the lattice spacing I. 

In the figure, we present a first, simple application of the method, the observation of the finite temperature cross-over 
from the discrete Tonks gas to the Mott phase in a system of impenetrable spinless bosons. This system for which 
U = +00 was totally out of reach for our previous algorithm with Hartree states 0] but can be easily treated by the 
stochastic Gutzwiller algorithm by simply setting iV max = 1. We introduce the name 'discrete Tonks gas' to refer to a 
phase of impenetrable bosons on a lattice with a filling factor much smaller than unity; as the mean distance between 
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bosons is then much larger than the lattice spacing, we expect this phase to have large scale properties very close to 
the ones of the continuous model, hence the name. The Mott phase corresponds to the opposite limit of a unit filling 
factor. 



VI. SEVERAL SPIN STATES 



All the discussion of the previous sections can be easily extended to the case in which the particles can be in s different 
spin states. The Bose-Hubbard Hamiltonian is immediately generalized as: 

H = Hu + Hj = E % 1 * t ( r > ff ')*(r, ^) + E E T S ^' ¥ ^ f7 )*( r '' ^ ( 45 ) 

r a.a'—l r,r' a—1 

\l/(r, a) being the destruction operators for a particle at the site r in the spin state a. The matrix U a , a > quantifies the 
on-site interactions between atoms respectively in the a and a 1 states. The hopping is assumed not to affect the spin 
degrees of freedom and is quantified by the parameters J a . 

The straightforward generalization of the Gutzwiller ansatz (j2J to the case of several spin states has the form: 

\G:4>)=H E 0(r,n)H ^ ' |0), (46) 



r ni ...n s =0 cr=l v °' 

where the Gutzwiller wavefunction <p(r, n) now depends on the spatial coordinate r and on the s-components vector 
n = (n%, . . . , n s ) formed by the occupation numbers of each spin state. 

As in the spinless case, both the on-site interactions and the mean-field contribution to the hopping term are taken into 
account by the deterministic evolution of 0(r, n). The fluctuations around mean-field are described by the stochastic 
noise term. 

For the case of a real-time evolution, the deterministic term reads: 
dt r ^ XJ aa > n a (n a > - 5 a ^,) 



+E Y 6 '" 1 ( * (r>) ^*( r ' ■>-<*) + i>"('',") N^7TT*(r, n + e,))-( jjj S »'w*"^ »)«(<■', »)) ») 

a.r' ii. u' .a 

(47) 

All the sums over spin states run over a = 1 . . . s and e CT is defined as the s-component unit vector whose <r-th 
component is I and all the others are vanishing. The mean- field ip(r, a) is defined as: 



The noise term can be written as: 

d<Mr,n)=E(- *iy ) 1 [^k,.e tkr (v^0(r,n-e ff )-^(r, ( T)0(r,n)) + 

k,CJ 

+d^ CT e- lkr (VH7TT0(r,n + e CT )-^(r,a)0(r,n))] , (49) 
in terms of the independent and zero- mean complex random variables dB^ a such that: 



dB^dBv,,,, = (50) 



dB^dBl.y = 2d T 5 kM ,S^ a ,. (51) 
w(k, a) describes the particle dispersion in the lattice in the absence of interactions: 

d 

fiw(k, a) — J fj E cos(fciZj). (52) 

i=l 
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As in the previous section, these equations can be extended to the imaginary-time evolution simply by replacing 
% — ► — 4f . The finiteness of the statistical error can be proven in the same way as previously done for the spinlcss 
case. 



VII. CONCLUSIONS 

In this paper, we have introduced a reformulation of the bosonic many-body problem on a lattice in terms of an 
stochastic Gutzwiller ansatz and we have proved that the statistical dispersion remains finite at all times in the 
presence of a cut-off on the number of particles per lattice site. This reformulation allows to perform efficient Monte 
Carlo simulations for strongly correlated systems like the Mott phase or the Tonks gas. This new method can also be 
used to study hard sphere bosons in arbitrary spatial dimensions, a situation which was totally out of reach for our 
previous stochastic Hartree ansatz. 
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APPENDIX: ANALYTICAL CALCULATION OF THE CORRELATION FUNCTIONS OF 

IMPENETRABLE BOSONS 

The correlation functions of a gas of impenetrable bosons in one dimension can be computed by mapping the bosonic 
system onto a fermionic system. The bosonic operators (x) on the lattice are written in terms of fermionic operators 
^ f (x) by means of the Jordan- Wigner transformation pll : 

x—l 

i>{x) = exp Utx ^f( x ')^f( x '))^f(x). (A.l) 

x'=0 

In this way, the Hamiltonian (JTJ for the impenetrable (U — +00) Bose gas reduces to a free-fermion Hamiltonian of 
the form: 

H = - [ *f(^ + t)*ir(x) + (-if-^U^riL -l)+ h.c] . (A.2) 

x=Q 

Given the specific ordering given to sites in l|A.l|) . it is in fact necessary to isolate in (|A.2|) the term corresponding 
to the tunneling between the last and the first site of the lattice, whose amplitude changes sign depending on the 
parity of the total number of particles present in the system. This peculiarity of the periodic boundary conditions is 
well-known in the frame of exactly solvable one-dimensional models, see e.g. |23j . 

For any value of the total number of particles, the Hamiltonian (|A.2fl is quadratic in the Fermi field operators and 
describes a homogeneous system of non-interacting fermions. The single particle states are plane waves of wavevector 
k and their energy is given by: 

hw{k) = Jcos(M), (A.3) 

but the quantization condition on k depends on the parity of N. For odd values of N, the Hamiltonian l|A.2|) reduces 
to the usual one for a one-dimensional lattice with tunneling constant J, so that periodic boundary conditions on the 
single-particle wavefunctions apply and the allowed values for k are given by: 

k = "-j-n (A.4) 

Li 
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in terms of the integer < n < J\f. On the other hand, for even values of N, anti-periodic boundary conditions are 
to be applied, so that k is now given by: 



k= T {n+-). 



(A.5) 



As the total number of particles N commutes with the Hamiltonian, we can write the Grand-Canonical density 
operator of our system as the sum of the contributions of respectively the odd and the even values of TV: 



P Podd ~t~ Peven; 



where: 



and 



Podd 



1 r 



1 r 
2 



-0H e . 



-0H odd +i%N 



(A.6) 



(A.7) 



(A. 



which involve modified Grand-Canonical density operators with ftp replaced by (3p + in. 
The odd- A" Hamiltonian i? dd reads: 



L-21 



H dd = - [ $U X + 1)*f(x) + & F (0)$ F (L - I) + h, 
while the even-A^ one reads: 



L-22 



ffcven = j [ £ ^f( X + 0*F(l) ~ *^(0)* F (i - i) + h. 



x=0 



-UN, 



fiN. 



(A.9) 



(A.10) 



In order to obtain the expectation value (O) of an arbitrary observable O, we have to calculate quantities such as Tr[cr] 
and (0) a = Tr[Ocr]/Tr[(r] where a is one of the four density operators appearing on the right-hand side of I) A. 7(1 and 
(|A.8p . As the a are Gaussian operators of the form exp(A) with A being a quadratic (but non- necessarily hermitian) 
operator in the Fermi field operators ^f(x), one has: 



(All) 



with the cifc's being the single particle eigenvalues of A. Depending whether one considers the first or the second 
term in l|A.7(l . I|A.8(I . the eigenvalues at are equal to —f3(fku(k) — /i), or to —j3(hw(k) — /x) + «7r, with the values of k 
corresponding to the relevant quantization conditions l|A.4|) or (|A.5|1 . 

The calculation of (0} a = Tr[Oer]/Tr[er] has to be performed using the Jordan- Wigner mapping (|A.1(I to write the 
observable in terms of the fermionic operators and then applying Wick's theorem, whose validity can be proven to 
extend to the present case of a non-hermitian quadratic Hamiltonian. In particular, the first order correlation function 
g^(x) in one of the density operators a can be shown |22,|24[ to be given by: 



g w (x?0) = (&(0)9{x)) a 



det 



2(& F (y)4> F (y'+l)) a -6 



(A.12) 



y,y'—0...x — l 



Note that we have added a global factor 1/2 which was missing in |22|. 
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